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Abstract 

We present results from 2D radiation-hydrodynamical simulations of fully compress- 
ible convection for the surface layers of A-type stars with the ANTARES code. Spec- 
troscopic indicators for photospheric convective velocity fields show a maximum of 
velocities near Tcff ~ 8000 K. In that range the largest values are measured for the 
subgroup of Am stars. Thus far, no prognostic model, neither theoretical nor numeri- 
cal, is able to exactly reproduce the line profiles of sharp line A and Am stars in that 
temperature range. In general, the helium abundance of A stars is not known from 
observations. Hence, we have considered two extreme cases for our simulations: a solar 
helium abundance as an upper limit and zero helium abundance as a lower limit. The 
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simulation for the helium free case is found to differ from the case with solar helium 
abundance by larger velocities, larger flow structures, and by a sign reversal of the 
flux of kinetic energy inside the hydrogen ionisation zone. Both simulations show ex- 
tended shock fronts emerging from the optical surface, as well as mixing far below the 
region of partial ionisation of hydrogen, and vertical oscillations emerging after initial 
perturbations have been damped. We discuss problems related to the rapid radiative 
cooling at the surface of A-type stars such as resolution and efficient relaxation. The 
present work is considered as a step towards a systematic study of convection in A- 
to F-type stars, encouraged by the new data becoming available for these objects from 
both asteroseismological missions and from high resolution spectroscopy. 



Introduction 

Research on A-type has benefitted enormously from the introduction of high resolution spec- 
trographs with R > 100,000 and from advancements in asteroseismology through telescope 
networks and satellite missions. They provide a new way of probing our theoretical under- 
standing of these objects and address questions such as: how are convection and pulsation 
linked to each other, in a physical parameter region no longer completely dominated by the 
strong envelope convection found in cool stars such as our Sun? How does convective mix- 
ing in these stars affect diffusion processes? What physical mechanism let chromospheres 
disappear in early A-type stars? How much does convection influence observed spectra and 
how does convection interfere with determining stellar parameters from spectroscopy? This 
selection of physic al questions demonstrates that A stars provide a uni que physical lab oratory 
( Landstreet , I2OO4 ) underlying the richness observed in their spectra ( Adelman . 2004 ). 



In the following we report on the first results from a project on numerical simulations of 
A-type (and F-type) stars, motivated by the progress made by asteroseismological missions 
such as CoRoT and MOST and challenged by results from high resolution spectroscopy. The 
simulations ultimately aim at improving our understanding of the dynamical processes of con- 
vection and its interaction with pulsation and, through mixing, on diffusion processes that 
shape the appearance of A-type stars. We first provide a background on how convection has 
observationally been identified in A-type stars. We then summarise the theoretical explana- 
tions given for the peculiar convective line broadening found in A stars and review previous 
attempts made to model the surface convection zones in these stars. In the subsequent 
section we discuss the setup of the numerical simulations we have performed for mid A-type 
stars followed by a discussion of our results and our conclusions as well as an outlook on 
further work. 



Observational Background 



The ex istence of convection in the surface layers of A stars was first predicted bv lSiedentopf 
( I933I ). Initially, observational evidence for such convection zones came from a non-zero 
microturbulence parameter ^. Common values for ^ were f ound in the range of 2 ... 4 km s~^ 
for both norma l A and Am-type main sequence stars (see lLandstreetj|2007 for a review and 
Adelmanll2004l for examples and further references). This readi ly exceeds the microturbulent 
velocities found for their cool star counterparts (see lCraxJligSSl for an extended discussion of 
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<^ in cool stars). To find direct evidence for convective velocity fields in the shapes of spectral 
lines of A stars is much more difficult. The average projected rotational velocities of normal 
A-type stars is found to be > 120 km s^^ and even for Am stars t;oSinz values of less than 
10 . . . 15 km are the exception and no t the rule (cf. the tables and references shown in 



Adelman and also the recent study ofl Royer et al.ll2007l ). Hence, first evidence that the 



bisectors of spectral lines in early type (sharp lin e) stars do not resemble those observed for 
cool stars was found f or late F-type supergiants (iGray &d Toner , 19861 ) and the early F-type 
bright giant Canopus (*D ravins & Lindl . Il984 : Dravind . 198 7^). Stars in this parameter range 
were shown to feature line profiles with curved bisectors caused by extended depressions in the 
blue wing of spectral lines. This results in an asymmetric line shape inverted in comparison 
with line profiles of cool supergiants (and those of our Sun and other cool main sequence 
stars). First evidence for this feature to be present in main sequence A stars was reported by 
Grav iS^ Nage] (|l989[ ) for the case of HR178 (H P 3883, an Am star). 



This was followed by a study of ^L andstreet (1998) who found two Am stars, HR4750 
(HD 108642) and 32Aqr (HD 209625), to show large asymmetries and a depressed blue wing 
in profiles of strong spectral lines. Weak lines were found to have much smaller profile 
asymmetries or even none at all. The same effects were found to be slightly weaker in the 
case of the hotter A-type star HR3383 (H D 72660). A hot A-star with similar properties, 
HR6470 A (HD 157486 A), was identified bv lSilaj et all (|2005l ). The current status of obser- 
vational evi dence of velocity fields caused by surface convection in main sequence A stars is 
analysed in iLandstreet et all ([2009), who added another five stars to this sample. The newly 
identified objects include a mid A-type star with stron g line profile asymmetries as previously 
found only for Am stars tGravl . [l989t ILandstreet . 199S) . Based on the collected data, which 
involved also reobserving objects during the same night, in consecutive nights, in different 
runs, and in different wavelength ranges, they concluded that the properties of the line pro- 
files of these stars are robust. Furthermore, they concluded that the only explanation for the 
line profiles which is consistent with the observational data is the presenc e of photo spheric 
velocity fields not related to pulsation thus confirming the suggestion by ILandstreet (19981) . 
These velocity fields are present in both normal A and in Am stars, though the latter contain 
the most extreme cases, in agreement with their higher values of ^. It is interesting to note 
that chromospheric activity indicators were identified for A-type stars with T^s ^ 8200 K 
( Simon et al. . 20021 ) which roughly coincides with the effective temperature of those stars 
for which the largest profile asym metries and overall convective line broadening have been 



identified (jLandstreet et al.l. l2009l ) 



Theoretical Explanations and Previous Simulations 

The most straightfor ward explanation for the "inverted" or "reversed" line bisectors has 
already been given bv lCrav &l Tonerl ([l986) for the case of supergiants who suggested that 
this could be caused by hot colu mns of rising gas ou tshining the cold downflows. This was 
further developed by I Gray ( 19891 ) and iDravind ( 1990l ) who noted that the hot upflows need 



to cover only a small fraction to reprod uce such lin e prof iles. This idea was also discussed 
for the case of main sequence A stars bv lLandstreetl (|l998. ). 

Could this be recovered from non-local models of convection or radiation-hydrodynamical 
simulations? A first study of the coupled convection zones caused by ionization of H i (and 
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He i) as well of He ii was presented in Latour et al. ( 198ll ). who used a modal expansion of 
an anelastic approximation to the hydrodynamical equations to describe convective motions 
in a mid A- type star. They noted the two convection zones to be coupled by overshooting. 
Due to the large He abundance (Y=0.354) assumed in this work it is difficult to compare it 
with more recent calculations (the large value of Y might partially explain why they found 
the con vective flux to be larger in the He ii ionisation zone, a property not recovered in later 
work). IXiond (jl990l) presented solutons obtained from his non-local model of convection 
for an F-type, an A-type, and two late B-type stars. Two separate convection zones were 
found coupled by overshooting for the mid A-type star with a convective flux being larger 
for the phot ospheric convection zon e . A s olar He abundance (Y=0.28) was assumed, as in 
th e study of Kupka &i. Montg omerx] (12002), who u sed the fully non-local convection model 
of Canuto &i Dubovik ov (1998) and Canuto et al. ( 200lh to construct envelopes of A stars 
for an evolutionary sequence of a 2.1 M© star and for an effective temperature sequence at 
three different metallicities. Except for the hottest A stars the two convection zones were 
found to be coupled by overshooting with the photospheric zone dominating in mid and 
late A-type stars (in agreement with the numerical simulations by Frevtag| |l995l . see below). 
They also investigated photospheric overshooting and found a positive kinetic energy flux 
for the o bserv a ble la y ers of m i d and l ate A-type stars. Such a distribution agrees with the 
model of iGravl (|l989[ l: lDravin3 (|l99Ql ): [L7ndstreej (|l998h . although no synthetic spectra were 
computed at that point. 

The first radiation-hydrodynamical (RHD) simulations of a mid A-type sta r in tw o spatial 
dimensions (2D) with 49 x 42 grid points were presented in ISofia &. ChanI (|l984'l (radial 
resp. vertical coordinate listed first here and in the following). However, an impenetrable 
upper boundary condition had to be placed at optical depths r of 0.22 to 2/3 in those 
simulations. This put the boundary layer right into the zone of partial ionisation of H i 
with its accompanying density inversion and pushed vertical velocities to zero in that layer. 
Consequently, the lower convection zone due to ionisation of He ii was found to dominate 
even at a solar He abundance and no predictions of photospheric velocities could be made. 
The s t ably s tratified layers in-between those two zones were found to be fully mixed. In 
Gigaj (|1989[ ) a 2D RHD simulation for the parameters of Vega was presented. The model 
had 30 X 36 points, but extended to much lower optical depths of logr —2.4. No He 
ionisation was included in the equation of state and simulation time had to be limited to one 
hour of stellar time. Vertical oscillations were found in the photosphere instead of a pattern 
of up- and downflows. Interestingly, line profiles cal culated using tha t simulation were found 
to be bent bluewards. However, a comparison with Landstreei (jl99& ) reveals that although 
the resulting bisector looks similar to those found for strong lines in HD 72660, a strong 
asymmetry i s also found fo r weak lines contrary to observations, as corroborated by the 
analysis of Landstreet et al.l ([2009), who included another four stars in that Teg range of 
8900 K to 10000 K in their analysis. Given the short simulation time and the reported 
periods of 5 to 15 min the observed oscillations could also have been transient. This is also 
implied by resu lts of the first extended study of A stars by means of 2D RHD simulations by 
Freytad ( 19951 ) (see also lFrevtag et al.l [l996). who found oscillations to be damped even for 
a Tcff ~ 9000 K and a \og{g) = 3.9, i.e., at the blue edge of the S Set instability stri p, while 
oscilla tions were clearly excited in an even cooler model. Most of the simulation in iFrevtag 
(|l995l) were carried out for a typical resolution of 65 x 62 grid points, with a few simulations 
doubling horizontal resolution, and one case reaching 95 x 182 points for a short run. A 
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solar chemical composition was assumed throughout and partial ionisation was accounted for 
as well. The case of non-grey radiative transfer was considered for one model and found to 
have little influence on the model structure and the flow itself. The convective (enthalpy) 
flux was always found to be larger in the H i driven convection zone compared to the deeper 
one driven by He ii (if included in the simulation box) for models with Toff < 9000 K with 
the two convection zones always connected by overshooting. A characteristic property of the 
simulations were nevertheless rather high velocities in the lower convection zone with vertical 
root mean square (wrms) maxima often around 1 km s^^. These maxima are related to strong 
vortices found in the lower convection zone. No synthetic spectra were computed from those 
simulations. 

A surprising change to these developments cam e finally with the f ir st RH P simulations 
in 3D for two mid A-type main sequence stars by iFrevtag &i StefFen ( 20041 ) (with a Tcff 
of 8000 K and 8500 K for a log(g) of 4.4). The grey approximation and a solar chemical 
composition were assumed. 90 to 110 grid points were used along the vertical direction 
and 180 points along the horizontal ones. Similar to the previous 2D RHD simulations and 
solar granulation simulations they followed a box-in-a-star approach in which a limited volume 
located near the surface of the star is considered for the calculations. Cont rary to most of the 



preceding 2D simulations the extent of the simulation boxes considered bv lFrevtag &i StefFen 



(j2004f l were chosen such as to contain at least several up- and downflow structures along 
any horizontal direction, which explains the larger number of grid points. Nevertheless, 
their simulations revealed the usual granulation pattern, with some subtle differences. In 
the end this lea d to line pr o files v ery much looking like solar ones in contradiction to the 
observations by Landstreetl ( 19981 ). The calculated line profiles were further analysed by 
StefFen et all (|2006l ). In a follow-up w ork a third 3D RHD sim ulation (for a Tes of 8000 K 



and a \og{g) of 4.0 ) was presented by iKochuk hov et al.l (12007 ) using again the CO^BOLD 
code as in .Frevtag & StefFenl (|2004l ) and StefFen et al.l (|2006l ). but this time with a slightly 
higher resolution (170 x 220^). In addition, part of the run was repeated using non-grey 
radiative transfer (based on binned opacities). Due to the latter the match for weak lines 
of stars such as HD 108642 (with actual stellar abundances used for the spectrum synthesis) 
improved. The same held for the cores of strong lines, but the wings of strong lines remained 
clearly discrepant both in shape, width, and bisector curvature. Only one more 3D RHD 
simulation of a cool A-star (Toff = 7300 K, \og{g) = 4.3, solar metallicit y, non-grey) with 
82 X 100^ grid points was presented in more detail bv lTrampedach ( 2004 ). though without 
computations of stellar spectra. However, in that simulation the H i and He ii zones are still 
connected to a single layer in a configuration described to be highly unstable. Further work 
in the literature has dealt with simulating interacting convection zones assuming idealised 
microphysics, but these results cannot be directly compared to stellar observations without 
invoking additional assumptions on the surface radiative cooling and the effects of chemical 
composition. 

The discrepancy found by Freytag &i StefFenl (2004) is surprising, since RHD simulations 
using the sa me code had been used to compute realistic solar spectra that match observed line 
profiles (see StefFen|[2007l and references therein). Obviously, the simulation results are in con- 
tradiction wi th the model of h o t an d narrow, rapidly rising columns of ga s. The observations 
discu ssed by iLandstreet et al" (|2009 ) indicate that chemical peculiarity ( Freytag &i StefFen . 



20041 ) can play only a limited role, as the line profile anomaly is found in both Am and normal 



A stars. This still cannot exclude He abundance to play a role (jKupkal . l2005l ). at least at the 
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level of explaining the differences between stars showing the most asymmetric and broadened 
profiles (of Am-type) and their counterparts am ong normal A stars. Effects of rotation might 
have to be consid ered, too (su ggested by Arlt in Frevtag &i Steffen| [2004. some consequences 
were analysed in iKupka 2005 ). This suggestion was also motivated by the frequent b inary 
nature of the ob s erved sharp line stars. However, the line profiles in Landstreet (|l998l ) and 
Landstreet et al. ( 20091 ) do not re semble at all t he fla t -bottomed line p rofiles of rapid ro- 



tators seen pole-on such as Vega ( Gulliver et al.l. 199 4: Hill et a 1. 1 , |2004| ). Since there is no 
indication from the data on the sample of Landstreet et a I., (,2009 ) that the binaries are very 
close or interacting, nor that there is any dependence of the profiles on the rotation period 
other than the ability of detecting them, rotation a ppears to be a less likely candidate to 
resolve this discrepancy. Rather, numerical resolution ( Frevtag &l Steffenl . 12004.) or the entire 
scheme used for the radia tive transfe r (Landstreet, priv. comm. 2004) may play a much more 
important role (see also lKupka 2005h . 



Numerical Simulations of A-type Stars with ANTARES 

To investigate the cause of these discrepancies between observation and numerical simula- 
tions we have decided to perfo rm RHP simulations of A-type stars with the ANTARES code 



( Muthsam et al. 2007, Mu thsam et al.ll20 09). This work is part of a more extended research 



project on convection in A- and F-type stars intended to study convection-pulsation inter- 
action and the properties of convection under physical conditions that lead to very different 
efficiencies for the transport of heat and for the mixing of fluid, respectively. 

As a first step we decided to have a closer look at the effects of helium abundance and of 
resolution, since both have not been investigated yet in sufficient detail. In addition, we also 
wanted to clarify the n ature of the extremely restrictive time st eps due to the short radiative 
cooling time scale ti-ad (|Frevtag &. Steffenl . l20o3 ISteffenl. I2OO7I ) (see Eq. (O). As ANTARES 



currently solves the fully compressible equations of hydrodynamics coupled to the stationary 
radiative transfer equation with a fully explicit time integration method, the ensuing time 
step constraints due to t^ad are severe and the resulting computational costs very high (10 to 
100 times higher than for a solar granulation simulation of comparable resolution according 
to [frevtag ii Steffen 2004). We have thus begun our study with RHD simulations in 2D to 
single out the most promising cases for 3D simulations and possibly suggest improvements 
to the time integration to reduce the computational costs of 3D simulations. 
With ANTARES we solve the RHD equations for fully compressible flows, 

|^-HV-[pu]=0, (1) 

V- [puu+p/] =pg + V-r, (2) 
V ■ [u(e + p)] = p(g • u) V • (u • r) + 0,.ad, (3) 



dpu 



de 
'Ft 

which are the continuity, Navier-Stokes, and total energy equation describing the conservation 
of mass density p, momentum density pu, and total energy density e. The later is the sum of 
internal (eint) and kinetic energy density, respectively, and u is the flow velocity. An equation 
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of state (EOS) relates eint and p to the temperature T and the gas pressure p (or rather, the 
sum of gas pressure and isotropic radiation pressure Prad> the meaning of p as used in the 
following). The EOS is taken from the most up-to-date OPAL tables (Rog ers et al. . 1996( ) 



and used also to compute all other required thermodynamic derivatives. We note that / 
is simply the identity matrix, r is the viscous stress tensor, and t and x denote time and 
space. The simulations discussed below are box-in-a-star simulations for a limited volume 
ranging from the upper photosphere to the stably stratified layers around T ^ 60,000 K to 
70,000 K, thus containing the surface convection zones inside the box. Since this region is 
small with respect to the stellar radius, a Cartesian (box) geometry is taken for the simulation 
box and the gravity g is a constant, downwards pointing vector equalling the surface gravity 
g. Horizontal boundary conditions are hence periodic and the geometry is plane parallel. The 
horizontal domain is chosen wide enough to fit several up- and downflow structures into the 
box at photospheric depth. The vertical boundary conditions are assumed to be impenetrable 
and stress-free (free-slip) with a constant energy flux assumed at the bottom of the simulation 
zone (and directly related to the Tcff defined for the simulation). The lower boundary is placed 
sufficiently far below the convection zones to have only limited influence (at most through 
reflecting waves). The upper boundary condition is placed at logTioss ^ —5 for similar 
reasons. The appropriate domain sizes are estimated from scaling the linear dimensions of 



solar granulation simulations by the pressure scale height I 
when comparing to the Sun implies 

1/Iq = Tcsp,QgQ/{Tcs.^Qpg) 



Hp - P/ipg) ^ T/Uig), which 



(4) 



for a star characterised by (Toff, ^,17), where p is the mean molecular weight. In optically 
thick layers the radiative (heating and) cooling rate Qrad = V • Frad can be computed 
from the diffusion approximation Fj-ad = — -K'radVT, where A'^ad is the radiative conductivity 
obtained from equation of state quantities and the Rosse land opacity Kross- For the latter, we 
use tabula ted OPAL opac i ties |lglesias &d Rogers! . Il996l ) extended at low temperatures with 
tables by iFerguson et al.l (|2005V both computed for the solar mixture of Greve sse & Noelsl 
(|l993ll . "or the photospheric layers a radiative transfer equation must be solved for accurate 
computation of Qrad- This computation is extended into deeper layers to allow a smooth 
transitio n to the diffusion app roximatio n. Details on this p rocedure as used in ANTARES are 
given in Obertscheider (2007) and in , Muthsam et al.l (|2009.) . We note that even in the grey 



approximation the results of this approach are different from the plain diffusion approximation. 
In ANT ARES the stationary limit of the radiative transfer (int ensity) equation is assumed (as 
in iFrevtaa .iggsl iFrevtag &i SteffenI I200I iKochukhov et al.l [2007.) and the resulting Q^ad 
is computed from an integration of the intensity obtained for a finite set of rays from this 
procedure. The grey case differs from the non-grey one by just integrating along each ray 
for a single frequency band rather than for multiple ones, but contrary to the diffusion 
approximation this accounts for inhomogeneity of cooling and heating in a convective stellar 
photosphere. We also point out here that the assumption of stationarity (instantaneous 
radiative transfer) leaves p, pu, and e as the dynamical variables of the system which have 
to be followed by the time integration during the simulation. 

To perform simulations which are close to an observational case we have chosen Tcff = 
8200 K and a log(.g) = 4.0 a s basic ph ysical parameter s, as they are representativ e of 
HD 108642 (|Kupka et al.l 12004 note that lLandstreeill998l and iLandstreet et al.ll2009l sug- 
gested a Tcff 100 K lower and a \og{g) 0.1 dex higher than those values, but such differences 
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are well within the observational uncertainties). Since diffusion has been suggested to sub- 
stantially change th e helium content of the upper en velope during the main sequence life 
time of A-type stars (jVauciair et al.l . ri974HRichard et al.. 2001; Michaudl l2004l ). we consider 
two extreme cases of helium abundance: a solar one and no helium at all. A solar metallicity 
mass fraction of Z ^ 0.02 was taken in both cases (this value is a compromise, because the 
underabundance of light elements relative to their solar values lowers the absolute value of 
Z while the opacity resulting from the different chemical composition contributing to Z is 
increased in the photosphere of Am-type stars; see K upka et al.ll2004l for further details on 
the abundances of HD 108642). The hydrogen mass fraction X is scaled up from 0.7 to 0.98 
in the case where helium is absent {Y = 0). 

For the simulations presented here we have chosen a coarse grid with a single refinement 
zone embedded inside the domain covered by the coarse grid. Since the simulations are in 2D, 
the horizontal boundary conditions of the coarse grid are periodic and because the fine grid 
is used to improve resolution near the surface, it has the same horizontal extent as the coarse 
one. Hence, it has periodic boundary conditions as well. Vertically, the fine grid is located 
in the interior of the coarse grid which provides the vertical boundary conditions for the fine 
grid. The coarse grid itself has the closed vertical boundary conditions already described 
above. For the case with solar helium abundance, the coarse grid has 150 x 200 points which 
span a domain of 32.2 x 40.3 Mm^ and yield a resolution of 216 x 201 km^. The fine grid has 
205 X 200 points and is located vertically between about 1.1 Mm and 12.1 Mm as measured 
from the top of the simulation box. This provides a resolution of 54 x 201 km^, i.e., a vertical 
refinement by a factor of 4. For the plots shown below the vertical axis is shifted by —3.8 Mm 
to have the downwards pointing vertical axis be labelled with zero where the horizontally and 
temporally averaged temperature equals Toff. For the simulation without helium, a coarse 
grid of 130 x 210 points is used spanning a domain of 26.1 x 41.7 Mm^ and thus providing 
a mesh size of 202 x 199 km^. The fine grid is located vertically between about 1.8 Mm and 
13.1 Mm as measured from the top and has 225 x 210 grid points. Hence, the resolution on 
the fine grid is 51 x 199 km^, which again implies a vertical refinement by a factor of 4. For 
this simulation the vertical axis is shifted in the plots shown below by —4.6 Mm. As a result, 
in plots showing both simulations the vertical axis value of zero indicates the average location 
where T = Toff, i e., the optical surface, for both cases. For the radiative transfer, we have 
used 12 rays and a single bin for Rosseland mean opacities, i.e., the grey approximation. 

Apart from setting up the simulations and some test runs which were done on single and 
two CPU core configurations most of the simulation runs to evolve both the solar helium 
and the helium-free case were done on 8 CPU cores using the MPI capabilities of ANTARES. 
The code was slightly modified to allow grid refinement zones to cover the entire horizontal 
extent of the model, as the fine g rid was previously use d to have higher resolution only 
in a specific domain of interest (cf. iMuthsam et al. 20071). These enhanced pos sibilities to 



use fine grids are now a standard feature of the code ( Muthsam et a I. . 20091 ). Another 



improvement added to the code and motivated by simulations of A-type stars was its ability 
to work with arbitrary chemical compositions (and thus combinations of X ,Y, and Z), since 
the opacity and equation of state tables used have some restrictions on the combinations 
of these parameters. The improved microphysics interface is now also a standard module 
of the code. After initializing the simulations, all production runs were performed on the 
P0WER5 p575 of the RZG in Garching. The initial conditions for the s imulations wer e taken 



from ID models for HD 108642 calculated with the LLmodels code (jShulyak et al.l. l2004l ) 
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using the convect i on mod el bv lCanuto &l Mazzitelli (1991) and opacity distribution function 
tables by 'Kurucz' ('1993a') and 'K uruczl (,1993bi) . For the ID models ^ = 4 km : 



assumed (Landstreet, 1998; Kup ka et a I 



was 

■2004l a nd T^g and log(.g) were the same as in the 
simulations (values taken from IKupka et a! .1120041. The flow in the simulations is started by 
a lack of perfect hydrostatic equilibrium of the stratification introduced by the differences in 
the equation of state and opacities used, the change in grid spacing between the ID model 
and the simulation, and the different numerical methods used in both codes. In addition, we 
add a small perturbation in the mass density field by creating a random distribution for the 
density perturbations in Fourier space. The values of this smooth noise function are added at 
each point to the average for the horizontal layer as obtained from the ID model. The noise is 
equally distributed over all modes (wave numbers), but its magnitude is scaled as a function 
of depth to produce the largest density fluctuations where the ID convection model predicts 
the largest temperature fluctuations. This avoids the introduction of long-lived perturbations 
in the radiative region near the bottom of the simulation box, which would cause unnecessarily 
long relaxation times (we return to this topic further below). 



Results 

in the following we discuss the main results of our 2D RHD simulations for both the case 
with solar helium abundance and for the case with zero helium abundance. Some results 
from experimental runs with different initial conditions and resolution but otherwise identical 
parameters are included as well. 



Relaxation and oscillations 

Both simulations were run for several sound-crossing times until convective motions had 
sufficiently developed, which is indicated by a growth in vorticity and in total kinetic energy. 
For the simulation with solar He abundance notable growth of convective motions occurred 
after about three sound-crossing times (time for a sound wave to cross the entire simulation 
box vertically, here tgc = 1371.6 s). After t/tgc 5 convective motions dominated the 
contributions to kinetic energy and experienced a phase of exponential growth. Figure [1] 
shows the velocity of the centre of gravity of the entire simulation volume as a function 
of time. The oscillations triggered by the stratification being out of perfect hydrostatic 
equilibrium are prominent and are almost purely vertical since the horizontal component is 
totally negligible. There is only weak intrinsic damping. We attribute this to the quality of the 
advection scheme implemented in ANTARES. We removed these oscillations by introducing 
a damping term at t/tgc ^ 6.1 for the vertical vel ocity component, as is routinely done so 



during relaxation of solar granulation simulations ( Trampedachl . Il997l ). Due to the short 
time scale chosen the oscillations are rapidly damped which allows us to turn off this artificial 
damping after t/tgc ~ 7.9. At this point the vertical oscillations contribute < 3% to the total 
kinetic energy of the flow and shortly afterwards the total kinetic energy starts dropping for 
the first time after its rapid phase of growth and at least the upper convection zone itself is 
essentially relaxed. The snapshots in Figs. [5l [7] and [8] have hence been taken at t/tgc 8. 
The equivalent holds for the averages in Fig. [3] and [H 

An alternative way of visualizing vertical oscillations is to look at the motion of the 
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Figure 1: Vertical (solid line) and horizontal (dashes) velocities of the centre of gravity 
for the simulation with solar helium abundance as a function of time. The horizontal 
component is magnified by a factor of 10"^. On the top x-axis, time is normalised by the 
sound-crossing time. Arrows indicate the beginning and the end of the artificial damping 
phase. 



photosphere and define the velocity Wph as the horizontally averaged vertical motion of the 
layer for which T = Tcs (the photospheric radius) at a given point in time. For both 
simulations this velocity is about an order of magnitude larger than that of the centre of 
gravity. Figure [2] shows the time development of Vph for the simulation with zero helium 
abundance. Intrinsic damping is just slightly more efficient than in the previous case. Artificial 
damping of vertical motions is applied from t/tsc ^ 4.4 to t/t^c ^ 13.9 onwards (here, 
^sc — 1150.9 s). As we have used a smaller damping rate for this second simulation the 
damping phase is longer and smoother than for the previous case. The development of 
convection also takes longer: strong convective motions set in only after t/tgc ^ 7 and they 
dominate over the vertical oscillations in terms of kinetic energy after t/tgc ^ 9. While the 
convective motions increase in strength, the contributions by oscillations drop until they reach 
a level of < 3% at t/tgc ~ 14. Once the flow can evolve undamped, some oscillatory motions 
seem to rise again. This is more obvious for Wph than for the centre of gravity, despite t'ph is 
also much more influenced by events taking place at the surface such as the formation of shock 
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fronts, which contribute to the non-sinusoidal shape of the function visible in Fig. [2] With 
kinetic energy increasing no longer exponentially the simulation slowly reaches equilibrium 
at t/tsc ~ 18. The snapshot in Fig. [5] and the averages in Fig. [4] have been computed for 

t/tsc > 17. 
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Figure 2: Temporal evolution of photospheric velocity for the simulation with zero helium 
abundance. On the top x-axis time is normalised by the sound-crossing time. Arrows 
indicate the beginning and the end of the artificial damping phase. 



We notice that a proper choice of the initial perturbations is more important for sim- 
ulations of surface convection in A stars than for solar granulation simulations. An earlier 
model for the case of solar helium abundance, which had been perturbed with an additional, 
sinusoidal velocity field instead of a random perturbation in density, showed strong resiliency 
to "forget" this pattern in the zone of He ii ionisation even after several sound-crossing 
times. That simulation run was hence discarded. This behaviour is quite different from the 
solar case where the same kind of perturbation is ra pidly removed by the convective flow, for 
instance, in the simulations of lObertscheider (|2QQ7h . 
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depth [Mm] 

Figure 3: Temporal and horizontal average of actual (V) and adiabatic (Vad) temperature 
gradient for the simulation with solar helium abundance. The regions of partial ionisation 
of H I, He I, and He ii are indicated. Layers unstable according to the Schwarzschild 
criterion are indicated in grey. 



Structure of the convection zone 



In Figure[3]we show the dimensionless temperature gradient V and the dimensionless adiabatic 
gradient Vad for the simulation with solar helium abundance. The superadiabatic peak in 
the zone of partial h ydrogen ionisation is a bout 3.5 times steeper than in a simulation of 
solar granulation (cf. [Rosenthal et al. 19991 for the latter). The upper convection zone is 
about 2.52 Mm deep as obtained from the definition of the Schwarzschild criterion applied 
to the horizontally and time averaged gradients. The upper zone is mostly driven by the 
partial ionisation of H i. It is extended due to the partial ionisation of He i which takes 
place close enough to unite both regions into a single convection zone. With 5.7 Mm the 
second convection zone, which is caused by partial ionisation of He ii, is much larger in 
vertical extent. This is mostly due to the increase of the local pressure scale height. The 
resulting mean structure is very similar to that one found in one dimensi onal models of stellar 
atmospheres and s t ellar envelopes fo r that part of the HR diagram (cf. Vauclair et al.lll974 , 
Latour et al.lll98l[ lLandstreelj|l998l ). because the high radiative losses of convection keep 
the stratification of the entire envelope in these stars close to that one of a purely radiative 
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Figure 4: Convective (enthalpy) and kinetic energy flux in units of the total flux, shown 
as solid and dashed line, respectively, for the simulation with solar helium abundance 
(black thin) and without helium (grey thick). The kinetic energy flux is multiplied by a 
factor of 10 for better visibility. 



mode l (this has also been foun d by lFreyta j 1995 . Freytag et al. 19961 Kupka & Montgomeryl 
2002llFrevtag &l Steffen|[2004l ). 

By comparison, the simulation with zero helium abundance has just a single convection 
zone due to the partial ionisation of H i. Its extent into the photosphere up from where 
on average T = Toff is marginally smaller than in the previous case (0.23 Mm instead of 
0.25 Mm). With a total extent of 1.56 Mm when following the same definition as used 
above, it is more shallow than its counterpart containing helium, in spite of a larger local 
scale height due to a lower /i, Eq. (j?]). This is a direct result of the lack of He i ionisation 
extending the zone deeper inside (note the behaviour of Vad in Fig. [3]). However, with a 
maximum of 1.6 the gradient V itself is flatter and the resulting superadiabatic gradient is 
only about 2.3 times steeper than in solar granulation simulations. 

Differences are also visible in the flux distribution. Figure |4] displays the enthalpy (con- 
vective) flux and the flux of kinetic energy for both simulations in units of the total 
(input) flux after subtracting the contribution due to the mean flow (i.e., due to the re- 
maining vertical oscillations). That contribution is negligible for Tk, but can contribute a 
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Figure 5: Snapshot of temperature fluctuations ([T'] = K) relative to their horizontal 
mean for the simulation with solar helium abundance. Extreme values have been chopped 
for easier visualization. The two black dashed lines indicate the region where grid refine- 
ment has been applied. The Tioss = 1 'isosurface' is denoted by a solid line. Streamlines 
indicate the direction and magnitude of the local flow field (the mean vertical velocity 
has been subtracted to improve visibility of the weakest vortices). 



relative flux of up to 0.5% to Fi^, if the averaging time is short, i.e., t < tsc- The higher 
enthalpy flux found in the simulation with no helium is hardly significant since it could be 
the signature either of a more relaxed simulation or of an effectively stronger convection due 
to a lower mean molecular weight. The same could hold for the higher velocities found for 
the helium free case, though we consider them more likely to be caused indeed by a lower 
/X. Nevertheless, and interestingly enough, we note that in the simulation with no helium 
i^k has the opposite sign for the upper half of the convection zone. This feature remains 
robust even when averaging for t/tgc ^ 3 and also when averaging over one third of that 
time scale for different subintervals. This inversion in the sign of i^k, which does not appear 
in simulations of solar convection, occurs in our simulation with the most shallow convective 
zone. In this case with a really thin convective layer heating from below and cooling from 
above act at the same time within a pressure scale height, while in a thicker zone, such as 
in the model with solar helium abundance, heating and cooling are more disconnected (cf. 
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Figure 6: Snapshot of temperature fluctuations ([T'] = K) relative to their horizontal 
mean for the simulation without helium. Scaling of T', domain boundaries of the grid 
refinement zone, the location of the optical surface defined at Ti-oss = 1, and the local 
velocity field (with mean vertical velocity subtracted) are indicated as in Fig. [51 



Moeng &i Rotunnd [l990l for a discussion of heating from below and cooling from above in 
a meteorological context). We als o note that at the ch osen scale th e fluxes in the He ii 
ionisation zone are not visible (cf. StefFen et al.l 20061 and IStefFenI 20071 ). These small fluxes 
require a more careful inspection of the simulation. 

Indeed, if we look at Fig. [5] we see that the strongest vortices created by the flow are 
located in the upper convection zone or just slightly underneath it and thus are close to the 
region of the strongest convective driving. The large downdraft at a horizontal coordinate of 
23 Mm manages to penetrate just a little bit into the lower convection zone, but at that time 
no strong vortices have developed in that region. We recall our previous observation that 
if vortices are seeded there as an initial condition, for instance by a sinusoidal perturbation, 
they remain there for a long time. This is both due to the weak convective driving (which 
cannot rapidly destroy a vortex or push it into a stably stratif ied layer) and due to the general 
longevity of vortex patches in 2D (cf. ^ Muthsam et al. 2007 ). Temperature fluctuations are 
found largest at the optical surface (the extrema have been truncated in the plot for a better 
contrast and can be four to five times larger while root mean square fluctuations are twice 
as large). The optical surface is corrugated mostly around rapidly evolving upflow regions. 

In Figure[6]we see quite a similar overall pattern for the simulation without helium. The 
horizontal scales are somewhat larger as expected from Eq. (|4]). Vortex patches are also 
strongest inside the convection zone itself, but a larger number of them is found further 
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below than in Fig. O We think that this is a result of the longer time evolution of this 
simulation. Eventually, also in Fig. [5] downflows from above will form new vortex patches 
in the lower convection zone or patches located in the stably stratified layer in between will 
reach it. Thus, if we run the simulation with solar helium abundance significantly longer, we 
may eventually see an increase in Fh. However, this may not be more realistic, as the vortex 
patches in 2 D just happen to assem ble near the bottom of a simulation of a strongly stratified 



medium (cf. iMuthsam et aTbOOTf l. Since 3D simulations started from 2D ones may inherit 



these properties for quite some (simulation) time, we think that future long-term simulations 
are the only way to accurately compute the enthalpy and kinetic energy fluxes in the lower 
convection zone of models for A stars with a helium content large enough to drive that zone. 

Shock fronts and their development 

As we can conclude from the small enthalpy fluxes found for both simulations shown in Fig.|4] 
co nvection i s rath er inefficient in t r anspo rt ing heat in mid A-type s tars, a r esult also reported 



by^'Freytag (l99^. 'Frevta g et aTl (Il996ll iKupka &i Montsomervl (|2002l ). iFrevtag &i StefFenI 



20041 .StefFen et al... (2006h and lSteffenl (|2007t l~ Consequently, the mean temperature gra 



dient stays close to the radiative one and hence, contrary to RHD simulations of solar granu- 
lation, a density inversion appears in the layer where H i is partially ionised ( Frevtaei . 119951 ). 
This is a consequence of the lower mean molecular weight which in turn is caused by the 
doubling of the number of particles contributed by ionised hydrogen relative to neutral one. 
In the Sun efficient convection and thus a much flatter V is counterbalancing this ef fect of 



partia l ionisation (cf . the simulations of M. Ste ffen and F.J. Robinson compared in iKupka 



I2OO9I: the model bv lCanuto &i Mazzitelli 1991 is an exception, since it predicts inefficient 



convection for the photospheric layers of the Sun which are located on top of a rapid transi- 
tion to efficient, adiabatic convection dominating the bulk part of the convective envelope). 
Figure [7] shows a snapshot for our simulation with solar helium abundance. The density 
inversion is prominent and the ensuing Rayleigh-Taylor instability enhances the convective 
instability caused by high opacity and a low Vad (see Fig. [3]for the latter). 

Note that the inversion layer is interrupted by a complex feature at a horizontal coordinate 
of 33 Mm. In this region the optical surface is more corrugated (solid line in the temperature 
plot in the middle panel) and supersonic velocities are found there, too (bottom panel). This 
feature is the result of two counterrotating vortices, which form at some time t/tsc ^ 7.5. 
At t/tsc ~ 8 they have created an extended region just underneath the photosphere with a 
temperature lower than its horizontal average. The low temperature fluid is supplied by the 
strong downdraft between the vortex pair which reaches supersonic velocities in its inflow 
area within the photosphere. At t/t^c ^ 8.2 hot gas from the neighbouring upflows has 
mostly covered the cold downflow region leaving a dense spot of cold material in its centre. 
As a result supersonic flow appears both in the two photospheric inflow areas and in the 
downdraft itself. From now on the structure begins to collapse emitting a series of shock 
fronts. At t/t^c ^ 8.4, which is shown in Fig. [71 there are no more supersonic inflow regions. 
Rather, supersonic velocities are found in a small front in the upper photosphere (bottom 
panel), visible also in T (band structure trailed by material of lower T underneath) and in T' 
in Fig. [5l which shows the state of the simulation 68 s before that one illustrated in Fig. [T] 
The streamlines in Fig. [5] also demonstrate that the structure is no longer dominated by a 
pair of roughly equally strong vortices. Supersonic flow is also found just underneath the 
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Figure 7: Mass density p, logarithm of temperature T ([T] = K), and local Mach number 
in the simulation with solar helium abundance prior to the formation of a large shock 
front. The isoline for which Tioss = 1 is denoted in the middle panel by a thick solid 
line. Isolincs for logTross = — 4,— 3,— 2,-1, +1 are denoted as dotted lines. Note the 
density inversion near the optical surface in the top panel. Regions of supersonic flow 
are surrounded by a bright solid isoline where Ma = 1. 



patch of dense gas extending into the photosphere (Fig. [7]). 

Once this gas patch has bumped onto the layers just below it, a large shock front is 
created which horizontally extends over one third of the simulation box and also spreads over 
almost one third of the photosphere contained in the simulation volume. This can be seen in 
Fig. [8l at ^ 0.25 tsc after the state shown in Fig. [71 The snapshot shows the moment when 
the remainders of the previous front, which is shown still moving upwards in Fig. [7] and is 
reflected by the upper boundary condition shortly afterwards, are absorbed by the new front. 
This explains the extended high Ma number region preceding the new front (we recall here 
that a shock front is a (near) discontinuity in the hydrodynamic variables, which does not 
necessarily have to move at supersonic velocities). As the upwards moving material has much 
higher density, it immediately absorbs the downwards moving one with very little changes to 
its physical state. Note that the shock front is optically thin, although with a steep gradient 
in optical depth (isolines of logTross — —3 and —4 nearly merging with each other). At the 
same time the anomaly in the density stratification has nearly disappeared. Comparing with 
Fig.[7]we can see that the density in this region has dropped dramatically (in fact to average 
horizontal values) leaving a low temperature region which is reheated by the upwards moving 



18 Effects of resolution and lielium abundance in A star surface convection simulations 




10 20 30 40 



[Mm] 

Figure 8: Mass density p, logarithm of temperature T {[T] — K), and local Mach number 
in the simulation with solar helium abundance 343 s (0.25 tsc) after the time step displayed 
in Fig. [71 Note the large shock front in the photosphere spreading more than 1/3 of this 
layer within the simulation box. It is characterised by large jump in temperature and a 
relative increase in the optical depth by more than an a factor of 10. Isolines have the 
same meaning as in Fig. [T] 



shock front. We note that similar events also appear in the simulation run without helium and 
smaller events appear more frequently than the extreme case illustrated by Fig.[7]and Fig. [8] 
These violent events have no counterpart in simulations of solar convection performed with 
the ANTARES code (Muthsam et al.l. l2007l ). Even though they are intermittent in nature, 
they are common enough (cf. the front visible in Fig. [6] in the upper right corner of the 
simulation volume) to be important statistically. They may even contribute to the shape of 
ine p rofiles and provide an efFicient source of heating up a chromosphere (cf. .Simon et aT] 



2002 



Resolution and computation of the radiative cooling rate Qrad 

Already for the case of solar granulation the most difFi cult quantity to resolve in equations 
(IT])-([3|) is the radiative heating and cooling rate Qiad ( Nordlund &i Stein , 199X1 ). It directly 
influences the large, energy carrying scales by driving the cooling (and reheating) of gas in 
the photosphere. Quantities directly related to Qiad such as the radiative flux Fiad and the 
superadiabatic temperature gradient V — Vad or the profiles of T and p are readily resolved 
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at vertical grid spacings of 12 . . .2 5 km, but to obtai n an e qually smoothly resolved Qrad 
requires a 10 times higher resolution ( Nordlund &i Steinl . Il99li) . even though T, p, and Fiad 
appear reasonably well converged at a resolution of ~ 12 . . . 25 km. Apparently, at that 
resolution the simulation runs benefit from some averaging effects. We also recall that in 
regions where the diffusion approximation holds, we have Qrad = ~div(i4'i.adVT) which 
involves a second derivative of T. This explains why resolution requirements for Qiad are 
higher than for mean structure quantities [p, T, ...) and their first derivatives (Fiad, V, 
...). 

How well do we resolve Qiad in A-type stars? Analysing our simulations for the coarse grid, 
which at 200 km resolution is equivalent to a solar granulation simulation with a resolution 
of ~ 50 km, the answer is: at that grid spacing not at all. Near the surface Qrad often has a 
triple-peak structure covered by about 15 grid points. At the same time on the fine grid there 
is only a double peak-structure ranging just 5 such grid points, i.e. 20 points on the fine grid, 
whereas the maxima of these peaks have less than 1/3 of the size found on the coarse grid. 
Physically, the double peak-structure originates from the radiative cooling layer at the bottom 
of the photosphere (which is also present in RHD simulations of solar surface convection) and 
a heating layer just underneath it. The latter is created by the partial ionisation of H i and 
the inefficiency of convective transport in the superadiabatic layer coinciding with the entire 
surface convection zone (see Fig. [3]). Once the simulation contains a fully developed surface 
convection zone, the detailed structure of the peak becomes more complex, as heating and 
cooling layers are no longer confined to a thin zone (cf. Fig.[5]to Fig. [8]). Since the coarse grid 
solution is discarded by ANTARES at each time step for a domain where a fine grid solution 
is available, the wrong solution on the coarse grid has no impact on the development of the 
fine grid. But it nevertheless implies a warning concerning simulations of mid A-type stars 
with lower effective resolution for Qrad than the one presented by our fine grid solutions. We 
consider even our current simulations on the boarder line of resolving the overall structure 
and magnitude of Qiad, because each of the very sharp peaks is represented by only about 
5 points. By comparison, at a vertical resolution of 16 km (corresponding to about 60 km 
for our present case) a solar granulation simulation would have the main (negative) peak 
covered by 15 points which resolve Qiad smoothly everywhere except near its extremum 
( Obertscheider . lioOTi ). Thus, A stars require a 3 to 4 times higher effective resolution at the 
stellar photosphere to properly resolve Qiad- This is also in agreement with the fact that V 
is steeper than in the Sun by just about that factor (Fig. [3]). 

There is a second problem associated with Qiad- Th e radiative time scales in the pho- 
tospheres of A star s are much shorter than in the Sun (jFreytagI, 1 19951 : Freytag & Steffenl 



2004t ISteffenl . 120071 ). since the opacity at their surface is lower while temperatures are higher. 



For optically thick layers the lower densities near their surface play a role, too. Hence, the 
time sca l e for relaxing a temperature perturbation of arbitrary optical thickness by radiation 
(|Spiegell.[l957h . 

m;^['-T '''''''' T J ' 

becomes smaller than the hydrodynamical time scales in a simulation. This includes the time 
scale of sound waves travelling a grid distance h, i.e., isound — h/cs. In Eq. ([5]), Cv is the 
specific heat at constant volume, k is the opacity, a is the Stefan-Boltzmann constant, and 
a perturbation of size I with k = 2tt/1 is assumed. For the optically thick case t^d converges 
to the time scale of radiative diffusion, tdis = 3{Kp/k)^Cv/{lQKaT^) ^ P/x< but remains 
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larger than zero for the optically thin case defined by Kp/k ^ 1 (using the specific heat 
at constant pressure, Cp, in the definitions does not change the argument). Comparing the 
time scales igound and tj-^d for the case h ~ I reveals the shortest radiative relaxation time 
scales in the problem. It is important to note that although for an A-type star the sound 
speed Cs is larger than for the Sun in the layer where T — Tcs, trad become s even smaller . 
For the Sun at usual resolutions of 20 .. .30 km for granulation simulations ( Kupka . 20091 ) 



tsound ^ ^lad, while for an A-type sta r with an equivalent resolution trad ~ 0-0 1 . . .O.ltgound 
for layers around the optical surface (jFrevta 3, 11995I: iFrevtag &. Steffenl . 12004 ^. 



As a consequence, RHD simulations of mid A-type stars, which use a p urely explicit time 



ntegra tion metho d, are lim i ted to time steps At < trad- This was noticed in lFrevtag &i StefFen 



(|2004l ) and also in Frevtag (1995). Since ANTARES currently uses a purely explicit time in- 
tegration scheme for Eqs. (IT])-([3]), it is subject to the same restrictions. The most severe 
limitations implied by ([5]) are found for optical depths 1 ;$ t < 10. There, tdiff is already a 
useful approximation of trad, whence At < h^/x as for the heat equation. Consequently, a 
grid refinement by a factor of 4 in that region requires a 16 times smaller time step. This 
explains the extremely small At of 5 ms and 3.4 ms for our 2D simulations with solar helium 
abundance and without helium, respectively. We note that the coarse grid time s t eps o f 



about 0.08 s and 0.05 s are comparable to the < 0.2 s reported in lFrevtag &i StefFenI (j2004l ) 
— these differences in At show that the 2D simulations presented here require computational 
efforts just slightly smaller than state-of-the-art 3D simulations at lower resolution h. But 
are such small At unavoidable, if Qrad is computed on the same mesh as the hydrodynamical 
variables? To check this we have computed the evolution time scale of the independent 
hydrodynamical variables p and e. If we know both variables at the grid at a time step n, 
i.e., p*^"^ and e("\ for any integration method we require that At < Cp'^"^(9p'"-'/c't)^^ 
and At < Ce(")(9e("V^^)"^ where C is a constant less than 1, typically 0.1, to be able 
to predict the new state of the system at the time step n + 1 (otherwise, even with fully 
implicit methods, iterative solvers may not converge, instabilities can occur, etc.). If we take 
C — 0.1, this is equivalent to requiring that at each grid cell the solution should not change 
by more than 10% during an integration step. Inspecting the simulation at three subsequent 
time steps we can easily estimate the time derivative and evaluate these inequalities. We have 
done this for the first time steps of the simulation and for a number of time steps spread over 
the entire duration of the run for solar helium abundance. To interpret the results we have 
also evaluated the convective flow time scale tc — h/ max((u^)^/^) as well as tg = h/ maxcs, 
and finally At < Ce'"^ (qJ;"^)^^ to see, if Qrad operates on a shorter time scale than any of 
the mechanical terms in Eq. ([3]). 

In Table[T]we show the results for the third time step, for a time step after initial relaxation 
through radiation, and for a time step during the generation of shock fronts described in the 
previous subsection. Time steps between the second and the third example show a rather 
continuous transition between these two. We note that only during the first time steps the 
radiative heating and cooling (e/Qrad) constrains the temporal evolution by providing the 
shortest time scale (even, if C = D). However, already after a few seconds the temperature 
perturbations have been smoothed out to an extent that sound speed and its associated time 
scale ts set the restrictions for the time evolution of the dynamical variables of the system. In 
the last time step shown supersonic velocities in the photosphere finally have led to tc < tg. 
We conclude that except for the first few hundred simulation time steps totalling just a few 
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Table 1: Time scales of the simulations with solar helium abundance. Grid 1 denotes 
the entire coarse grid, grid 2 the fine grid, grid 3 the coarse grid without the domain for 
which the fine grid is available. C = 0.1 and D = 0.25 in all calculations. 



time scale 


grid 1 


grid 2 


grid 3 




t ^ 0.06 


s 






1.27 


0.51 


1.27 


Dt, 


18673.50 


2204.89 


36166.30 


6 mm(p/pt) 


1.48 


31.11 


1.48 


C min(e/et) 


0.56 


0.12 


1.49 


Cmin(e/(3i.ad) 


0.40 


0.12 


17.97 




t = 3.62 


s 




DU 


1.27 


0.51 


1.27 


Dt, 


371.56 


23.12 


587.73 


Cmin{p/pt) 


25.60 


7.60 


25.60 


Cmin(e/et) 


7.33 


2.91 


20.33 


Cmin(e/(3rad) 


0.84 


2.15 


64.36 


t = 


11442.92 s ( 


~ 8.3 tsc) 




Dt, 


1.27 


0.51 


1.27 


Dt, 


1.86 


0.49 


8.05 


Cmin{p/pt) 


0.80 


0.53 


6.30 


Cmin(e/et) 


0.81 


0.53 


7.35 


Cmin(e/(3rad) 


0.40 


1.12 


13.59 



seconds, which are subject to the assumed initial conditions and random perturbations applied 
to the latter, the temporal evolution of the system is governed by the hydrodynamical time 
scales. This holds even for Qrad itself, which implies changes on the total energy on a similar 
time scale as hydrodynamical processes (note that the overestimation of Qiad on the coarse 
grid leads to a slightly shorter time scale, but this part of the solution is discarded by using 
the fine grid solution — outside the fine grid domain the coarse grid solution imposes no 
restrictions). 

Thus, mathematically Qrad is just a stiff term in a differential equation making the whole 
problem at least in principle solvable on the evolution timescale of the dynamical variables 
of the system by a properly designed implicit integration method. The restrictions in At are 
solely due to high wave number components k contained in Qrad, which represent radiative 
transfer over one or a few grid cells. This transfer indeed occurs on short time scales trad, but 
does not govern the evolution of the system itself, as it takes the much longer time Ce/Qrad 
to substantially change the energy content of a grid cell radiatively. A natural approach is to 
consider an operator splitting technique and integrate only Qrad by an implicit method. This 
strategy is used also in simulations of convection i n rota ting spherical shells to model the 
lower part of the solar convection zone ( Clune et al.l . ll999l ). But here the difficulty is that the 
coefficients in Qrad are neither constant, nor linearly dependent on (p, e). Plain subcycling 
for the integration of Qrad (multiple radiative transfer steps per hydrodynamical time step) 



22 Effects of resolution and fielium abundance in A star surface convection simulations 



alone is inefficient, since Qrad as a whole does not cause rapid changes to e, but only its 
high k components evolve that way, while a consistent update of its coefficients is non-trivial. 
Thus, we consider it more promising to analyse higher order methods with proper damping of 
small but rapidly evolving components or filtering methods for their suitability to accelerate 
RHD simulations of stellar convection at high resolution. This would bring 3D simulations 
resolving the radiative heating and cooling at the surface of mid A-type main sequence stars 
from a supercomputing application to the realm of high performance depart ment computers , 
as ha ve been used in solar granulation simulations with ANTARES (see iMuthsam et al.ll2007 , 
2QQ9h . 



Conclusions and Outlook 

We have presented 2D RHD simulations performed with the ANTARES code for a mid A-type 
star for two extreme cases of helium abundance, a solar one and a helium free composition. 
The simulations differ from previous works by a higher resolution and the application of 
high (5*'') order advection schemes running stably for these simulations without the need 
to introduce artificial diffusion (see Muthsam et al.i r2009 for further details). The quality 



of the advection scheme was demonstrated by the need to introduce artificial damping of 
vertical velocities over several sound-crossing time scales of the simulations to remove vertical 
oscillations introduced by the initial conditions. After that phase oscillations driven by the 
flow itself are found to reappear. The mean structure of the convection zone is close to 
a radiative one, as found in previous RHD simulations in 2D and 3D. The most interesting 
differences found between the two cases we have considered include an inversion of the sign of 
the flux of kinetic energy as well as higher velocities and larger flow structures, each of them 
observed for the case with zero helium abundance. This can be understood in terms of the 
smaller extent of the convection zone due to the absence of partial ionisation of helium and 
the lower mean molecular weight of a helium free mixture. Since the evolution time scales 
of the surface convection zone are sufficiently short, we expect these results to be robust. 
As a note of caution we point out the lack of strong vortices found for the zone of He ii 
ionisation for t he c ase with solar abundance, which is in contrast to the results reported in 
Frevtad (1995) and Frevtag et a 1. 1 (l996. ). From further simulation runs with different initial 



conditions we have found that at least for the lower (He ii ionisation driven) convection zone 
the simulations should be performed over very long times to ensure they no longer depend on 
the initial perturbations (vortices seeded initially may just remain while they may take a long 
time to grow on their own). As soon as the surface convection zone is sufficiently developed 
large shock fronts are emitted into the photosphere. These result from the interaction of 
vortex pairs with the density inversion near the surface. We expect that a similar mechanism 
could work around strong downflow regions in 3D simulations as well. The large shock fronts 
might affect spectral line profiles and provide a mechanism of heating chromospheres in mid 
A-type stars. Finally, we have analysed the importance of resolution to properly compute the 
radiative heating and cooling rate Q,ad and have shown that the time step restrictions on a 
hydrodynamical solver implied by computing this quantity on a fine mesh are just those of a 
classical stiff term in a differential equation. Properly designed implicit integration methods 
should thus be able to accelerate RHD simulations of this class of stars. 

This is most important since the two simulations shown here required 8.3 million time 
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steps for the fine grid solutions resulting in a total of ^ 13, 000 CPU core hours (on P0WER5 
1.9 GHz CPUs). A 3D RHD simulation at the same resolution h is hence a supercomputing 
project (we recall that this is a consequence of the At < h'^/x dependence of an explicit 
solver — a larger h as used in previous works would make the simulation more affordable, 
but this should be traded in only, if Qrad can be computed sufficiently accurately). Such a 
project is feasible with the ANTARES code which has been successfully run with good scaling 
on up to 1024 CPU cores at the P0WER6 machine of the RZG in Garching. However, this 
would still restrict us to very few simulation runs. We thus think it is important to work on 
the implementation of proper integration methods, since this would provide benefits also for 
simulations of other types of stars once high enough resolution is demanded. Moreover, it 
would bring 3D RHD simulations of mid A-type stars with a resolution h as presented here 
into the realm of high performance department and university computers, which is already 
the case for equivalent simulations of F-type to M-type main sequence stars. 

We are working on a version of ANTARES with open vertical boundary conditions which 
avoids the reflection of waves and shock fronts. This will make the simulations more suitable 
for the study of convection-pulsation interactions and increase the stability of the code, since 
it avoids extreme conditions occurring when a front hits the upper boundary. We intend to use 
long-term simulations with this new version of ANTARES to compute synthetic spectra for 
comparisons with observations, study convection-pulsation interactions, and probe non-local 
models of convection, followed by 3D simulations as a final step. 
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